% default is 1 second simulation of RS neuron.
% For a chattering neuron try:
%p.C=50;p.vr=-60;p.vt=-40;p.k=1.5;p.a=.03;p.b=1;p.c=-40;p.d=150;p.vpeak=25
%izhikevich_neuron(500,p)
function spikes=izhikevich_neuron(stim_current,params)

if ~exist('params')
    C=100;vr=-60;vt=-40;k=3;
    a=0.01;b=5;c=-50;d=100;
    vpeak=35;
else
    C=params.C;
    vr=params.vr;
    vt=params.vt;
    k=params.k;
    a=params.a;
    b=params.b;
    c=params.c;
    d=params.d;
    vpeak=params.vpeak;
end
T=1000;tau=1;


spikes = 0;
n=round(T/tau);
v=vr*ones(1,n); u=0*v;
I=[stim_current*ones(1,n)];

for i = 1:n-1
    v(i+1)=v(i)+tau*(k*(v(i)-vr)*(v(i)-vt)-u(i)+I(i))/C;
    u(i+1)=u(i)+tau*a*(b*(v(i)-vr)-u(i));
    if v(i+1)>=vpeak
        v(i)=vpeak;
        v(i+1)=c;
        u(i+1)=u(i+1)+d;
        spikes = spikes + 1;
    end
end
plot(tau*(1:n),[v' u']);

